clear; if ispc; dbstop if error; end; rng('default');

input_folder = [fullfile(pwd, '..') filesep];

%% load data and demand, mc, and fc estimates
list_variables = {'para_nonlinear','cov_var','randomdraw','dim_price','price_coeff_opt','model','logincnorm','rc_option',...
            'yr','craftdummy','quantity','geomkt','store_county_id','loginc', 'est_corr', ...
            'owner_name2','fixed_brand','fixed_craft_brands','list_brand','brewer_name','brand_descr','pop_profit','prod_attributes',...
            'craft_mkt','list_county','price_coeff','setup','cutoff_month_count', 'N_boot'};


load([input_folder 'fc_covariates2016.mat'], list_variables{:});

load([input_folder 'gms_cs.mat']);

addpath([fullfile(pwd, '..', '1demand') filesep]);
addpath([fullfile(pwd, '..', '2mc') filesep]);


%% setup
cost_reduction_ratio = 0; %0 (no merger efficiency), 1 (merger efficiency)


cf_firm_type = {'entry', 'merged', 'incumbents, other'};

base_est_years = 2016;
cf_est_years = 2016;
change_prod = true;

check_sim_fc = true;
starting_point_option = [];

n_fc_para_sampled = 10; % # of fc parameters sampled 
n_fc_sim = 5; % # smln draws of fc shocks
fixedcost_option = 'conditional';


% market size cutoffs for defining small/median/large markets
pop_profit_cutoff = [0 1e5 5*1e5 Inf];

%% sample from the confidence set of the fixed-cost parameters
sample_from_cs;
if use_scope    
    para_in_state = -para_vec(1, :);    
    para_craft = -para_vec(2, :);   
    para_fc0 = -para_vec(3:5, :);
    para_sig = para_vec(6:8, :);
    para_first = -para_vec(9:11, :); 
    para_mkt_sig = zeros(1, n_fc_para_sampled);
elseif est_corr
    para_in_state = -para_vec(1, :);    
    para_craft = -para_vec(2, :);   
    para_fc0 = -para_vec(3:5, :);
    para_sig = para_vec(6:8, :);
    para_first = zeros(1, n_fc_para_sampled);
    para_mkt_sig = para_vec(9, :); 
else
    para_in_state = -para_vec(1, :);    
    para_craft = -para_vec(2, :);      
    para_fc0 = -para_vec(3:5, :);
    para_sig = para_vec(6:8, :);
    para_first = zeros(1, n_fc_para_sampled);
    para_mkt_sig = zeros(1, n_fc_para_sampled);
end

%%  preparation
%(1) assign products and indices by owner
list_owner = unique(owner_name2);
N_owner = length(list_owner);

prod_name_by_owner = cell(N_owner, 1);% products for each owner
prod_id_by_owner = cell(N_owner, 1);% index of these products in list_brand
brewer_name_by_owner = cell(N_owner, 1);% brewers for these products

fixed_prod_ind_id = find(ismember(list_brand, fixed_brand));

brewer_est_yr = brewer_name(yr==base_est_years);
brand_est_yr = brand_descr(yr==base_est_years);
owner_est_yr = owner_name2(yr==base_est_years);

for no = 1 : N_owner    
    brand_no = unique(brand_est_yr(ismember(brand_est_yr, list_brand) & ismember(owner_est_yr, list_owner(no))));
    prod_name_by_owner{no} = setdiff(brand_no, fixed_brand);    
    prod_id_by_owner{no} = find(ismember(list_brand, prod_name_by_owner{no}));    
    
    [~, brand_id_no] = intersect(brand_est_yr, prod_name_by_owner{no});    
    brewer_name_by_owner{no} = brewer_est_yr(brand_id_no);        
end
clear brewer_est_yr owner_est_yr brand_est_yr

%(2) select markets to simulate
[~, sortind] = sort(-pop_profit);
geomkt_select = craft_mkt(sortind);

%(3) define merger
% rank the total quantity sales of the craft breweries
craft_brewer_list = unique(brewer_name(craftdummy & yr==base_est_years));
q_craft = nan(length(craft_brewer_list), 1);
for k = 1 : length(craft_brewer_list)
    q_craft(k) = sum(quantity(yr==base_est_years & ismember(brewer_name, craft_brewer_list(k))));
end
[~, qsort_ind] = sort(-q_craft);
clear craftdummy


merge_owner = [setdiff(...
    unique(owner_name2(ismember(brewer_name, craft_brewer_list(qsort_ind(1:5))))), ...
    {'BOSTON BEER COMPANY', 'SIERRA NEVADA BREWING COMPANY', 'MILLERCOORS', 'HEINEKEN',  'AB INBEV'}); 'AB INBEV'];
divest_brands = [];
merged_craft = find(ismember(list_owner, setdiff(merge_owner, 'AB INBEV')));


if ~isempty(merge_owner)
    merge_owner_id = find(ismember(list_owner, merge_owner));
else
    merge_owner_id = [];
end

clear brewer_name

%(4) define firm move sequence
total_q = nan(N_owner, 1);        
for no = 1 : N_owner            
    total_q(no) = sum(quantity(ismember(owner_name2, list_owner{no})));
end
[~, owner_seq] = sort(-total_q);
clear quantity owner_name2 total_q included

%% organize data by markets for parallel computation in cf simulations
n_smln_mkt = length(geomkt_select);
ind_yr = ismember(yr, cf_est_years); clear yr

% define mktsize dummies
pop_profit_by_mkt = cell(n_smln_mkt, 1);
mktsize_by_mkt = cell(n_smln_mkt, 1);
county_order_by_mkt = cell(n_smln_mkt, 1);
prod_entry_by_mkt = cell(n_smln_mkt, 1);
fixed_prod_entry_by_mkt = cell(n_smln_mkt, 1);
sim_pcoeff_by_mkt  = cell(n_smln_mkt, 1);
loginc_by_mkt = cell(n_smln_mkt, 1);


for m = 1 : n_smln_mkt  % do not use parfor here, which requires a lot of memory
    pop_profit_by_mkt{m} = pop_profit(craft_mkt==geomkt_select(m));
    mktsize_by_mkt{m} = [pop_profit_by_mkt{m}<pop_profit_cutoff(2), ...
        (pop_profit_by_mkt{m}>pop_profit_cutoff(2) & pop_profit_by_mkt{m}<pop_profit_cutoff(3)), ...
        pop_profit_by_mkt{m}>pop_profit_cutoff(3)];       

    ind_m = (geomkt==geomkt_select(m) & ind_yr);  
    county_order_by_mkt{m} = find(ismember(list_county, store_county_id(ind_m)));  

    brand_tab_m = tabulate(brand_descr(ind_m));        
    count_brands_m = cell2mat(brand_tab_m(:, 2));        
    prod_entry_by_mkt{m} = brand_tab_m(count_brands_m>cutoff_month_count, 1);

    ind = ismember(list_brand, prod_entry_by_mkt{m});        
    fixed_prod_entry_by_mkt{m} = ind(fixed_prod_ind_id);

    sim_pcoeff_by_mkt{m} = unique(price_coeff(ind_m,:), 'rows');        
    loginc_by_mkt{m} = unique(loginc(ind_m, :), 'rows');                
end


clear brand_tab_m count_brands_m pop_profit geomkt list_county store_county_id brand_descr list_brand price_coeff loginc

%% simulate counterfactuals
% draw a base for fc shocks
% note: draw more than n_fc_sim to ensure that we have n_fc_sim fc draws
% after weeding out draws that cannot rationalize the observed product
% entry
prod_sim_fc_all = cell(N_owner, 1);
for no = 1 : N_owner
    prod_sim_fc_all{no} = rand(length(prod_name_by_owner{no}), n_fc_sim*20); 
end

disp('CF simulation start')


for np = 1 : n_fc_para_sampled
    prod_attributes_np = prod_attributes;

    % pick a set of draws of the demand supply shocks from their asymptotic distribution   
    % i.e., draw on the third dimension of demand_resid
    nb = randsample(N_boot, 1);
    for j = 1 : length(prod_attributes_np)
        if isfield(prod_attributes_np{j}, 'demand_resid')
            prod_attributes_np{j}.demand_resid = prod_attributes{j}.demand_resid_boot(:, :, nb);
            prod_attributes_np{j}.mc = prod_attributes{j}.mc_boot(:, :, nb);
            prod_attributes_np{j} = rmfield(prod_attributes_np{j}, {'demand_resid_boot','mc_boot'});
        end
    end
    
    parfor m = 1 : n_smln_mkt  
        fprintf(1, 'np = %1.0f, m = %1.0f\n', np, m);
    
        % get standard deviation parameter and normalize other parameters
        sig_m = max(1, abs(mktsize_by_mkt{m}*para_sig(:, np))); % cap the std at $1/mkt/product

        para_profit_n = 1/sig_m;
        
         
        para_in_state_n = para_in_state(:, np)/sig_m;        
        para_craft_n = para_craft(:, np)/sig_m;        
        para_fc0_n = (mktsize_by_mkt{m}*para_fc0(:, np))/sig_m;        
        if use_scope
            para_first_n = (mktsize_by_mkt{m}*para_first(:, np))/sig_m;
        else
            para_first_n = para_first(np)/sig_m;
        end        
        if est_corr
            para_mkt_sig_n = para_mkt_sig(:,np)/sig_m;
        else
            para_mkt_sig_n = para_mkt_sig(np)/sig_m;
        end
        
        % define merger efficiencies
        cost_reduction = cost_reduction_ratio*para_craft_n;    

        % define the initial configuration
        config_init_by_owner = cell(N_owner, 1);        
        for no = 1 : N_owner            
            config_init_by_owner{no} = ismember(prod_name_by_owner{no}, prod_entry_by_mkt{m});            
        end

        % define the order of movement
        included = false(N_owner, 1);
        for n = 1 : length(owner_seq)
            if ismember('entry', cf_firm_type)
                if (sum(config_init_by_owner{owner_seq(n)})==0) && (~ismember(owner_seq(n), merge_owner_id))
                    included(n) = true;
                end
            end    
            if ismember('merged', cf_firm_type)
                if ismember(owner_seq(n), merge_owner_id)
                    included(n) = true;
                end
            end    
            if ismember('incumbents, other', cf_firm_type)
                if (sum(config_init_by_owner{owner_seq(n)})>0) && (~ismember(owner_seq(n), merge_owner_id))
                    included(n) = true;
                end
            end    
        end
        owner_seq_m = owner_seq(included);

        % define mkt-specific product attributes
        prod_attributes_m = prod_attributes_np;
        
        for j = 1 : length(prod_attributes_np)
            if isfield(prod_attributes_m{j}, 'demand_resid')
                prod_attributes_m{j}.demand_resid = prod_attributes_np{j}.demand_resid(geomkt_select(m), :);
                prod_attributes_m{j}.mc = prod_attributes_np{j}.mc(geomkt_select(m), :);
            end
            prod_attributes_m{j}.dist = prod_attributes_np{j}.dist(county_order_by_mkt{m});
         end
        
        % compute post-merger eqm

        [owner_profit1, owner_profit2, owner_profit, ...
            mean_p1, mean_p2, mean_p, ttq1, ttq2, ttq, ...
            cons1, cons2, cons, nprod1, nprod, is_same, prod_sim_fc, ...
            added_prod, dropped_prod, ind_cf]...
            = prod_equi_sim(config_init_by_owner, ...
            prod_id_by_owner, prod_attributes_m, brewer_name_by_owner, merge_owner, ...
            fixed_prod_ind_id, fixed_prod_entry_by_mkt{m}, change_prod, prod_sim_fc_all,[], n_fc_sim, pop_profit_by_mkt{m}, ...
            para_fc0_n, para_profit_n, para_in_state_n, para_first_n, para_craft_n, para_mkt_sig_n, cost_reduction, merged_craft, ...
            list_owner, sim_pcoeff_by_mkt{m}, setup, ...
            para_nonlinear, randomdraw, ...
            loginc_by_mkt{m}, logincnorm, starting_point_option, owner_seq_m, check_sim_fc, ...
            geomkt_select(m), np);
    end
end
